Elastic mechanics solution of thermal expansion of bi-material curved beam and its application to negative thermal expansion metamaterials

Thermal stress impacts various engineering fields significantly, such as aerospace and precision instruments. This adverse effect can be greatly reduced, if not eliminated, by the application of micro-thermal expansion materials, and bi-material beams are widely utilized in the design of micro-thermal expansion structures, thereby exhibiting great application potentials. The elasticity solution of bi-material curved beam under free thermal expansion has been proposed by scholars. Based on this solution, the simplified form is proposed in this paper, and extended to the case where the rotation angles at both ends of the circular arc are constrained under thermal loads. Besides, the geometric parameters and the nonlinear problems of the thermal expansion of bi-material curved beam are analyzed. In addition, a novel type of negative thermal expansion material has been designed by applying the bi-material curved beam to the tetra chiral and anti-tetra chiral materials. The proposed material has greater negative thermal expansion effect than the traditional tetra and anti-tetra chiral materials that are with straight beams.

www.nature.com/scientificreports/ counterclockwise from the x-axis of the curved beam. Based on the theory of elastic small deformation, Gonczi proposed the stress field and displacement field inside the curved beam, and the boundary condition when there is no bending moments at both ends of the curved beam 34 . For material i ( i = 1, 2): where A i , B i , C i , D 1 , D 2 , D 3 are the undetermined coefficients, σ ρi , σ ϕi represent the normal stress and the circumferential stress of material i, u ρi , u ϕi represent the normal displacement and the circumferential displacement of material i, respectively. ρ and ϕ represent the polar radius and the polar angle in the polar coordinates.E i , µ i , α i represent the modulus, Poisson's ratio and CTE of material i, respectively. T represents the change of temperature. It should be mentioned that the blue and red regions represent material 1 and 2, respectively. The boundary conditions are as follows: In Eq. (1), D 1 , D 2 , D 3 represent the rigid body displacements and rotations of a circular arc, and they are regarded as arbitrary constants without affecting the stress and strain field of the circular arc. Therefore, if the center of the circular arc is taken as the reference point, it can be set as D 1 = D 2 = D 3 = 0 . In this way, the number of undetermined coefficients in Eq. (1) is equal to the number of equations in Eq. (2), and the equations are solvable: where A 2 , B 2 , C 2 can be obtained by interchanging material parameters and a, c in A 1 , B 1 , C 1 .
By substituting Eq. (3) into Eq. (1), the displacement component at ρ = b can be obtained: Figure 2 shows the schematic diagram of the deformations at both ends of a curved beam. According to the geometric relationship, the relative deformation at both ends of the curved beam is: In this way, the equivalent CTE of the curved beam can be obtained:  (3) and remove the higher-order terms of h b , the following is obtained: By Substituting Eq. (9) into Eq. (7), the following is obtained: In this section, we will use the finite element method to compare the accuracy and scope of application of the derived simplified formula and the formula in Ref 31 [Eq. (12)].
The 8-node plane stress element CPS8R of the finite element software Abaqus was utilized. Take θ = 180 • as an example of numerical verification. Figure 3a shows the finite element model. The distance between the interface of the two materials and the center of the circle is b = 80 mm , and the thicknesses of the two materi- , as mentioned in previous studies 33 . The upper part represents Material 1 and the lower part represents Material 2, with a grid of 100 (circumferential) × 6 (radial) and geometric nonlinearity is not considered in this section. After testing, the calculation results converge under this density of grids.
The middle section of arc is fixed and a temperature load of 1 °C is applied to obtain the equivalent CTE. Figure 3b shows the deformation of the beam. Then, the accuracy of the obtained analytical solution is verified by checking the error of the analytical solution and the numerical simulation when h b changes.
With the thicknesses of two layers being equal ( h 1 = h 2 = h 2 and θ = 180 • ), the accuracy of the proposed simplified formula and Eq. (12) as h b varied is investigated. Figure 4a shows the equivalent CTE α calculated by the proposed simplified analytical method (Simplify, blue curve), Eq. (12) (Ref 31 , red curve), and the finite element method (FEM, scatters) when h b increases from 1/16 (i.e. h = 5 mm ) to 1/5 (i.e. h = 16 mm ). Further study revealed that the relative error of results calculated via Eq. (12) remains at around 5.5% while the results calculated via the proposed simplified analytical formula are highly consistent with numerical simulation results and the relative error remains at around 0, which indicates that the result is indeed the elasticity solution of the bi-material curved beam. It is worth noting that here, the relative error at a certain h b is defined as the absolute value of the difference between the results calculated via the analytical formula and the numerical simulation results divided by the largest absolute value of the FEM results. Therefore, the proposed simplified analytical formula can be utilized to obtain the equivalent CTEs of complex structures based on bi-material curved beams.  12) at h (bθ ) = 1 10 , that is, the slenderness ratio of the arc is 0.1. When θ is small (< 54°), numerical simulation results show that the α of the curved beam is positive, while the result of Eq. (12) is negative. Besides, for any θ , the relative error of results calculated via Eq. (12) remains at around 20% while the results calculated via the proposed simplified analytical formula are highly consistent with numerical simulation results and the relative error remains at around 0, which not only indicates that small θ is not the applicable scope of Eq. (12), but the accuracy of the proposed simplified formula is very high. This proves the advantage of the proposed simplified formula.

Discussions
Parameter analysis. First, the effect of the thickness of curved beam on its equivalent CTE is investigated. Figure 4a shows that the curved beam exhibits negative CTE, while NTE effect degrades as h b increases. This can be attributed to the increase of the flexural rigidity of the curved beam with the increase of its thickness. As h b increases, the curved beam is more difficult to bend, thereby its NTE effect is affected. Meanwhile, the absolute value of CTE ranges from 250-910 ppm/°C, which is an order of magnitude higher than that of component materials, specifically in the range of h b shown in Fig. 4a. Therefore, bi-material curved beams can not only achieve NTE but also magnify the sensitivity of thermal expansion.
Then, the effect of the central angle θ of curved beam on its equivalent CTE is investigated. Figure 4b shows the equivalent CTE of the curved beam decreases monotonously from positive values to negative values as θ increases from 40° to 180°. This can be attributed to the increase of the rotation angles at both ends of the curved beam during thermal deformation with the increase of radian, which leads to a larger NTE. Therefore, according to the analyses above, h b shall be minimized and θ shall be maximized to achieve the maximum NTE effect.
Also, the effects of the thickness ratio h 1 h 2 of the two materials on the equivalent CTEs of the curved beam were investigated. Figure 5 shows α versus h 1 h 2 at h b = 1 5 . Each curve corresponds to a different θ . As observed, α decreases then increases as h 1 h 2 increases. This phenomenon is due to fact that when there exists a large difference between the thicknesses of the two materials, the bi-material beam will be closer to the singlematerial beam which is PTE, and the NTE effect will not be good enough compared with the bi-material beam with a certain ratio of the thicknesses of the two materials. Therefore, to obtain the maximum NTE effect of bi-material curved beams, understanding how to obtain the minimum α is necessary.
Maximum NTE. In this section, the relationships between h 1 h 2 and other parameters are studied, in particular when the curved beam has the maximum NTE effect (minimum CTE). Take the derivative of Eq. (10) with respect to h 1 h 2 and let the expression equal 0, the analytical relationships between h 1 h 2 and other parameters when α reaches the extreme value can be obtained. When α reaches the extreme value, h 1 h 2 −2 is approximately proportional to E 1 E 2 and the proportional coefficient is not only dependent on h b and θ , but also independent from α 1 and α 2 , as shown in Fig. 6a, which indicates that the ratio of static moments of the two layers should be specified to achieve maximum NTE. Figure 6b shows the curves of a specific proportional coefficient corresponding to different h b and θ when maximum NTE effects are achieved, and these curves will become the guidance of the design of the maximum NTE of bi-material beams.  Nonlinear problems. The derivations of the above formulae are based on small deformation. For the finite deformation problem, the incremental theory can be used. Based on the updated Lagrangian description (ULD), the deformation of the first incremental step is obtained from the initial state, and then the deformation state of the first incremental step is used as a known condition to find the unknown quantities of the second incremental step, and so on until the final incremental step is done.
Assume that before the ith incremental step, the radius of curvature of the beam is b i−1 and the center angle is θ i−1 , then after the ith incremental step, the radius of curvature and the center angle of the beam become: The coefficients A 1 , B 1 , C 1 are shown in Eq. (9). It is important to note that these coefficients contain b, and need to be in incremental forms, too.
The average CTE of the arc is: It can be seen that the results obtained by Eq. (18) are in good agreement with the finite element results, and the maximum error is only 0.66%, which shows that the analytical formula in this paper is also applicable to the incremental theory. www.nature.com/scientificreports/ In addition, the incremental theory can also be applied to the range of material nonlinearity. Two kinds of materials with non-linear physical parameters, high chromium steel and austenitic stainless steel, whose average CTEs (compared with 20 °C) and elastic moduli are shown in Tables 1 and 2, are used for analyzing. Besides, the physical parameters corresponding to the temperatures that are not given in the tables can be obtained by linear interpolation. Considering that the CTEs given in the tables are relative to the initial configuration, that is, based on total lagrangian description (TLD), it needs to be transformed into the CTEs based on ULD, and the transformation formula is as follows: In the iterative process, the CTE and elastic modulus corresponding to the temperature of the final state are substituted into Eq. (18) for iterative calculation. It should be noted that in this algorithm, to obtain the equivalent CTEs of the temperature corresponding to a certain iterative step, it is necessary to take the temperature as the final state and use the corresponding physical parameters for another set of iterative calculations. Figure 9 shows the results of the comparison between Eq. (18) and the nonlinear finite element method. The geometric parameters of the beam and temperature increment step are b 0 = 80 mm , θ 0 = 90 • , h 1 = h 2 = 4 mm , T 0 = 20 • C , T= 0.5 • C . After the final incremental step, the radius of curvature of the beam is b final = 77.5 mm , and the central angle is θ final = 93.24 • . It can be seen that the results obtained by Eq. (18) are in good agreement with the finite element results, and the maximum error is only 0.8%, which shows again that the analytical formula in this paper is also applicable to the incremental theory.

Novel NTE metamaterials
Structure design. As shown in Fig. 10, two novel NTE metamaterials are designed in this paper by replacing the straight ligaments in the tetra chiral/anti-tetra chiral honeycomb with bi-material curved beams, where the blue parts represent Material 1 (PI) in the inner layers of the curved beams, and the red parts represent Material 2 (PMMA) in the ring and the outer layers of the curved beams. In the chiral honeycomb, the ligaments are divided into two segments that bend in opposite directions. For both metamaterials, when the ambient  www.nature.com/scientificreports/ temperature rises, the curvature of the ligaments will increase, making the rings closer to each other, which will approximately produce NTE effect on the whole. Besides, the bending of the bi-material beam during thermal deformation is the basis of realizing the NTE effect, which leads to the rotation of the rings in the chiral structure and the rotation amplifies the NTE effect of the whole structure. Therefore, the bi-material design makes more contribution to the NTE than chirality. In fact, for the chiral structure in Fig. 10b, when the circles shrink, the NTE effect will be weakened and when the circle shrinks to a point, the structure will degenerate into the structure in Fig. 7a, whose NTE effect is much weaker compared with the structure in Fig. 10b.
The equivalent CTE. First, the equivalent CTE of tetra chiral honeycomb is deduced. As can be seen in Fig. 11a, in the tetra chiral honeycomb, the radius of curvature of the ligament arc is b , the central angle is θ , and the radius of the ring is r.
According to the geometrical relationship, the chord length of circular arc is: and the distance between the two rings is: The change of the distance between the centers of the two rings after the temperature rises T is:  www.nature.com/scientificreports/ As can be known from "Derivation, simplification and verification of analytical formulas" section: Therefore, the equivalent CTE of the tetra chiral honeycomb is: Then, the equivalent CTE of the anti-tetra chiral honeycomb is derived. As can be seen in Fig. 11b, in the anti-tetra chiral honeycomb, the radius of curvature of the semi-ligament arc is b , the central angle is θ , and the radius of the ring is r . The difference between tetra and anti-tetra honeycomb lies in that the distance between the two rings in anti-tetra honeycomb is: The equivalent CTE of the anti-tetra chiral honeycomb can then be derived using the same method: Finite element model. The 8-node plane stress element CPS8R of the finite element software Abaqus is utilized to perform numerical simulation to verify the results derived in "Structure design" section. A single cell is used for analysis. The displacement loads in the y direction are applied to the midpoints of the cross section of the ends of the left and right ligaments of the structure, and the displacement loads in the x direction are applied to the midpoint of the ends of the upper and lower ligaments according to the antisymmetry. The temperature load of 1 °C is applied to the whole structure, and the equivalent CTE can then be obtained.
Validity and parameter analysis. In this section, the change of CTEs of chiral/anti-chiral honeycombs as the curvature of the ligament changes is studied when the radius r and distance L of the rings remain unchanged. Meanwhile, the validity of the results has been verified by the comparison with the finite element results.
As shown in Fig. 12, the analytical solutions are basically consistent with the finite element results, and the maximum relative error of chiral honeycomb does not exceed 9% while the maximum relative error of the Figure 11. The geometry parameters of the novel metamaterials: (a) tetra chiral structure; (b) anti-tetra chiral structure. www.nature.com/scientificreports/ anti-chiral honeycomb does not exceed 6%. With the increase of central angle, the NTE effects of both chiral and anti-chiral honeycombs increase, which means that the design method of replacing bi-material straight beams with curved beams can effectively increase the NTE effect of this kind of metamaterial. Besides, when the central angles are equal to 0, the structures will become the NTE structures proposed by Yu et al. 19 and it is obvious that the proposed metamaterials in this paper could obtain higher NTE values compared with those according to Fig. 12. The jagged shape of the curve for FEM in Fig. 12b is due to the great robustness at the junction between the ring and the ligament and specific optimization problems need to be further studied. It should be noted that AM still means the results from analytical methods and FEM means the results from finite element methods in Fig. 12.

Conclusion
In this study, the precise elasticity solution for the thermal expansion problem of bi-material curved beam is proposed based on the study in Ref 34 and the simplified result is proposed based on the assumption that the slenderness ratio is small. This solves the problem that the accuracy of the formula in Ref 31 at relatively large thicknesses (as compared to the radius) of the beam and small central angles are limited. The length of the simplified result is much shorter than that of the original result, and at the same time, it is still in great agreement with the numerical simulation result. Also, parameters and the nonlinear problems of the bi-material beam are analyzed based on the simplified formula during thermal deformation, and the conditions under which the bimaterial circular beam achieves maximum NTE are identified. Finally, a class of novel NTE metamaterials are designed by introducing bi-material curved beams into tetra/anti-tetra chiral honeycomb, and the equivalent CTEs of this kind of metamaterials are given by using the simplified formula proposed in this paper. The proposed analytical formula can effectively analyze the thermal expansion of bi-material circular beams and has great potential for applications.

Data availability
The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.